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Abstract We introduce a solvable model of randomly growing systems consisted of many independent 
subunits. Scaling relations and growth rate distributions in the limit of infinite subunits are analyzed 
theoretically Various types of scaling properties and distributions reported for growth rates of complex 
systems in wide fields can be derived from this basic physical model. Statistical data of growth rates 
for about 1 million business firms are analyzed as an example of randomly growing systems in the 
real- world. Not only scaling relations are consistent with the theoretical solution, the whole functional 
form of the growth rate distribution is fitted with a theoretical distribution having a power law tail. 
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1 Introduction 

Growth phenomena are generally highly irreversible dynamical processes far from thermal equilibrium 
40J. From the viewpoint of statistical physics it is an important new topic that growth rates of complex 
systems often show non-trivial similar statistical behaviors across fields of sciences. Fat-tailed distribu- 
tions of growth rates and a non-trivial shrink of variance as a function of the size are reported in various 
fields of sciences; business firms [UHMI-Sl] ' sa ^ es °^ P narmaceu tical products [ID] , circulation numbers 
of newspapers [23], population of migratory birds [IS], animals' metabolic rate fluctuations [17], the 
amount of scientific funding [23] , group size of religious activities [22 ! , population size of cities [T2J , 
country's whole economic activity observed by GDP |18j and the amount of exports and governmental 
debts |25j . Probability densities of growth rates in most of these examples are typically approximated 
by double exponential (Laplace) distributions or by power law distributions, quite interestingly not by 
a Gaussian distribution. 

Statistics on growth rates of business firms have a long history of study, and recently statistical 
physicists are involved in this topic. Gibrat postulated the "law of proportional effect" that the expected 
value of the growth rate of a business firm is proportional to the current size of the firm [L4"II2o] . The 
original Gibrat's assumption states that the variance of growth rates is independent of the size, however, 
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Fig. 1 Comparison of system size dependence of the probability density functions of growth rates for different 
values of a. a; The case of a = 1.5 (< ffj(i) 1 ' 5 >= 1) and b; a = 0.5 (< gj(t) ' 5 >= 1) in semi-log plots. In 
both cases the growth rate distribution of individual subunits follows a uniform distribution. N = 1 (black thin 
line), N = 2 (red dash-dotted line), N = 10 (blue broken line), N = 10 2 (bold green line), N = 10 s (purple 
dotted line), N = 10 4 (light blue dash-dotted line) and N = 10 6 (black broken line). 



data analyses of business firm activities show various types of variance-size relations. There are papers 
which support the original Gibrat's assumption [11(133] . on the other hand non-trivial fractional power 
laws are reported not only for business firms but also in many other phenomena |HIBU81 fl0"ll 1 2 II 1 5 l[T7l 
18,22,23,24,32,40 . Also, the country-dependence [21] and the transition from Gibrat's assumption to 
such power law decays are pointed out [37] . Various types of theoretical models of business firms have 
been introduced by physicists for better understanding of scaling properties of business firms from 
the standpoint of complex systems Pll3ll4l[7 l l2lI l [2 ^ [27 1 [3T 1 [5g i l4T | l4l] . however, there has been no unified 
theory which can explain all these basic properties simultaneously. 

In the next section we introduce a new type of random growth model of a complex system con- 
sisted of many independent subunits, and in section 3 we theoretically solve non-trivial scaling relations 
between the variance of growth rates and the number of subunits by introducing a kind of renormal- 
ization. In the limit of infinite number of subunits the distribution of growth rates is shown to converge 
to a stable distribution with a power law tail. The stable distribution and the generalized central limit 
theorem were established in mathematics about 80 years ago [91IT9]. however, these concepts were ap- 
plied mostly for theoretical models assuming scaling properties for phenomena such as turbulence [351 
136] . Validity of the theoretical results for growth rates is confirmed in section 4 by analyzing a huge 
database about business firms. This example may be the first real world application of an asymmetric 
stable distribution fitted for the whole scale range. The final section is devoted for discussion. 



2 The model 

We consider a system consisting of N subunits characterized by non-negative scalar quantities, {xj(t)}. 
For each subunit we assume the following random multiplicative time evolution [16) . which is known 
to be one of the basic process of producing power law fluctuations [2l?ll3"8] . (See Appendix A for brief 
review of random multiplicative process). 

x j (t + At)=g j (t)x j (t) + f j (t), (1) 

where gj(t) and fj(t) for j = 1, 2, • • • , N, are growth rates and external forces, respectively, both are 
assumed to be independent identically distributed random variables taking only positive values. In the 
case that the probability of occurrence of gj(t) > 1 is not zero and if < \og(gj(t)) >< 0, where < • • • > 
denotes the average, it is known that there exists a statistically steady state in which the cumulative 
distribution follows a power law [IS] . 

PfrxJccxJ", (2) 



Fig. 2 Confirmation of theoretical estimates of size dependence of the width of the growth rate distribution. 
The width of growth rate distributions are measured by repeating numerical simulations for various values of 
TV. The points are plotted in log-log scale for five cases: a = 0.5, 1.0, 1.5, 2.0 and 2.5, (black, red, green, blue 
and orange), simulation results(points) and theory (lines). For numerical simulations, the growth rate for each 
subunit, gj(t), distributes uniformly in the range (0, (a + 1) 1//q ] satisfying < gj ( t) a >= 1 (a = 0.5, 1.0 ,1.5, 2.0 
and 2.5). (a)Data shown are median of G 2 for numerical simulations and Eq. 10 for theoretical estimations. 
(b)Data shown are the interquartil e ra nge (IQR) of G for numerical simulations and the square roots of the 
asymptotic functional forms of Eq. [10] given in Table 1 for theoretical estimations. The theoretical curves are 
shifted along the vertical axis so that they overlap with the numerical simulations at TV = 10 6 in the panel (b). 
The IQR, which is the one of the most commonly-used robust measure of width of a PDF, is difference between 
the largest and smallest values in the middle 50% of a set of data. Fro m these figures, we can confirm that the 



medians of G(t; TV) 2 almost correspond to the width given by Eq. 10 and IQRs are proportion to asymptotic 
behaviour of this equation given by table 1. 



for large value of a;,, with the positive exponent determined exactly only by the statistics of the growth 
rate by solving the following equation, 

< 9j (t) a >=l. (3) 

Based on a renormalization point of view we pay attention to the sum of all subunits, X(t; TV) = 
2j=i x j{t)> which follows the same type of time evolution as that of the subunits, 

X(t + At; TV) = G{t; N)X(t; TV) + F(t; TV), (4) 



where F(t; TV) = fjW: an d the growth rate of the whole system is defined as 



It is easy to show that the mean value of growth rates is invariant, namely, < G(t; TV) >=< Qj(t) >= G. 

Properties of this model can be investigated by numerical simulation. Fig. [IJi shows an example 
of deformation of growth rate distributions for various values of TV in the case that Eq. ([3) is fulfilled 
with a = 1.5 observed in the statistically steady state realized for time steps larger than 10 6 . Here, 
the distribution of growth rates of subunits is given by an independent uniform distribution as shown 
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Fig. 3 Convergence of the normalized growth rate distributions in the case of a — 1.5. The growth rate for 
each subunit, gj(t), distributes uniformly in the range (0, (2.5) 2,/3 ] satisfying < Qjii) 1 ' 5 >= 1. The renormalized 
growth rate's probability density functions for systems with JV subunits are plotted for several values of JV. 
N = 1 (black thin line), JV = 2 (red dash-dotted line), JV = 10 (blue broken line), JV = 10 2 (bold green 
line), JV = 10 3 (purple dotted line), JV = 10 4 (light blue dash-dotted line), JV = 10 6 (black broken line), and 
the theoretical symmetric stable distribution, p(G; 1.5,0) (yellow dotted line). The inserted figure shows the 
cumulative distribution function of the positive part in log-log scale confirming convergence to the power law. 



by the case of JV = 1. In this figure, the additive noise term, fj(t), is set to be a positive constant for 
simplicity as we confirmed that the main results do not depend on the values of fj (t) except the case 
that it is identically 0. As known from this figure the distribution of growth rates of the aggregated 
system tends to shrink slowly to a delta function as JV goes to infinity. It is numerically confirmed that 
the same property of convergence to the delta-function holds for any distribution of growth rates of 
subunits if the growth rate distribution satisfies Eq. ^ with a > 1. 

Fig. [TJd shows an example in the case of a = 0.5. In this case the growth rate distribution stops 
shrinking for JV larger than 10 and it converges to a non-trivial distribution in the limit of JV goes to 
infinity. It is confirmed that this property is always observed if the value of a in Eq. ([3]) is between 
and 1. The distribution of growth rate in the limit of JV = oo depends on the functional form of the 
growth rate distribution for the subunits. 



3 Theoretical analyses 

We can theoretically evaluate the JV-dependence of the width of the PDF of G(t; n) by introducing an 
approximation of the random variables {xj(t)} which are known to follow a power law in the steady 
state. We introduce the following values as the measure of the width, G(t; JV) 2 . Here, we define G(t; N) 
as: 

G(i; JV) = G(t; JV) - G = ' J ■ (6) 

G(t; JV) 2 takes zero in the case that the PDF of G(t; N) is a delta function. Let it be a random variable 
following a uniform distribution in the interval of (0, 1]. Then, the distribution of the new variable, 
u —i/a^ f u ows a power law with exponent a. For uniform random variables, {u}, we can approximate 
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Table 1 Summary of generalized central limit theorems for growth rates. The value of a and the asymptotic 
functional form of Eq. llOl divided by a 2 , which corresponds to the width of the growth rate, for large system 
size, N, and the limit distributions. 7 = 0.577... is the Euler constant. 



N samples by a set of deterministic values, {1/N,2/N, ■ ■ ■ , 1}, so that the set of power law variables 
{xj(t)}, which follow Eq. can be approximated by the deterministic set: 

{JVV", (iV/2) 1 / Q , (iV/3) 1 /", • • • , (N/j) 1 / 01 , ■ ■ • , l}. Therefore, the typical sample of G 2 (t; N) is obtained, 



{EjliU/")- 1 '"}' 



Taking the average of G r (t;N) with respect to gj(t)(j — 1,2, • • ■ N) and applying the independence 
condition, < (g n (t) - G){g m (t) - G) >= a 2 S nm , we have 



< G r (t;N) 2 >= a 2 • J^ML ,s, 



where a 2 is the variance of the growth rates for the subunits and S nm is Kroneker's delta. Then, 

we can calculate the summations in Eq. (JsJ) , A rl ' Q EjLii _1 ' Q an d N 2 t a Y^f=i J _2 ^ Q j by applying an 
asymptotic expansion formula for the Riemann Zeta function, 

00 1 N 1 1 1 

C(A)=E -X=E ^+ (A _ 1) ^-i -2^ + — 0) 

Neglecting the third term and higher order terms in the right hand side of Eq. ([9| , we have the following 
approximation of Eq. ^ for < a: 

(G r (t,N) 2 )=a 2 - 2 ~ a ? . (10) 



These theoretical evaluations are checked numerically in Fig. [2j in which the widths of growth rate 
distributions for the aggregated system are plotted as functions of the number of subunits, N. The 
theoretical asymptotic functional forms in Table 1 fit quite well asymptotically for all cases. It should 
be noted that the ordinary standard deviation or variance is not a good measure for the width of the 
distribution of growth rates as these quantities arc affected easily by outliers. Instead of the standard 
deviation we introduce the median of G(t; N) 2 and the interquartile range (IQR) for characterization 
as the measure of the width of distribution, in Figs. [2|a) and (b) respectively. The IQR, which is the 
one of the most commonly-used robust measure of width of a PDF, is difference between the largest 
and smallest values in the middle 50% of a set of data. From Fig.^a) we can confirm that the medians 
of G(t;N) 2 almost correspond to the width given by Eq. 10 and from Fig.[2^b) IQRs are proportional 
to asymptotic behaviour of this equation given by table 1. 
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From Eq. ( 10 ) we find that the width of PDF of G(t; N) converges to following a power law of 
N in the case of 1 < a, while the width takes a finite value even in the limit of N —> oo in the case 
of < a < 1. At the marginal case of a = 1 the width decays to very slowly for increasing N. For 
2 < a the width decays inversely proportional to N, that is, the typical N dependence in the case of 
ordinary central limit theorem. These functional forms of N dependence are summarized in Table 1 in 
the second column. As known from this result Gibrat's assumption of constant variance is fulfilled in 
the case of < a < 1, and the non-trivial power law decays of width-size relations observed in many 
complex systems are realized in the case of 1 < a < 2. As shown in the first column of this table the 
range of a is characterized by the form of equality or inequality for the moment function of the growth 
rates, M(s) =< gj(t) s >. Derivation of these relations is summarized in Appendix B with the basic 
properties of the moment function. From this table we find that Gibrat's assumption holds for systems 
in which the mean growth rate is larger than 1, while the non-trivial power law decay of width-size 
relation for growth rate is expected in the situation that the mean growth rate of subunits is less than 
1. 

Next, we theoretically derive functional forms of the distribution of growth rates normalized by the 
width in the limit of N — > oo. We consider the following 3 cases according to the value of a. 

I: The case of < a < 1: The width of growth rate does not shrink to but it converges to a 
finite value even in the limit of N —¥ oo as known from Eq. ( |10| . This reason can be understood 
by a general property of power law distribution with exponent less than 1. In such a case the mean 
value < Xj(t) > diverges implying that some of samples in {xj(t)} take extraordinarily large values 
compared with others. So, both the denominator and numerator of Eq. ^ can be approximated by 



only finite numbers of extraordinarily large contributors, therefore, the value of Eq. ( 10 ) is finite even 
in the limit of large N. The distribution of growth rate shows the same property, namely, even in the 
limit of N — > oo, the limit distribution is determined only by a small number of large contributors, 
therefore, we cannot expect a universal functional form in this case. 

77. The case of 1 < a < 2: As indicated in Fig. 2 the width of distribution shrinks to in the limit of 
N — > oo we can expect existence of a universal limit distribution independent of the initial condition. 
In this case the average, < Xj(t) >, takes a finite value in the steady state, the denominator in Eq. (JgJ) 
can be roughly estimated as N < Xj(t) > for very large values of N. On the other hand, the numerator 
in Eq. @ is given by the sum of (</,•(<) — G)xj(t), in which ffj-(i) — G gives a coefficient taking either 
positive or negative sign randomly with respect to j, and Xj(t) follows a power law with the exponent 
a. Namely, the numerator becomes a summation of N independent identically distributed random 
variables that have both positive and negative power law tails with the exponent a. As the generalized 
central limit theorem can be applied to such a sum of random variables, the limit distribution of growth 
rate G , which is normalized by the width of the distribution around the mean value, is expected to 
converge to a stable distribution, which has the form of an inverse Fourier transform (see also 
Appendix C for a brief summary of the central limit theorem and the stable distributions) , 

p(G;a,/3)=^y exp [-ipG - \p\ a (1 - i^(p, a))} dp, (11) 

where the asymmetry parameter (3, which takes a value in the interval [—1, 1], and the function ijj(p, a) 
are given as follows, 



It is well known that the limit probability density, Eq. (11), has a power law tail with the exponent 
a just like the distribution of {xj(t)}. 

In Fig. [3] we confirm the validity of this theoretical result by a numerical simulation for the case of 
a = 1.5 and (3 — 0. The normalized growth rates for the system consisted of TV-subunits, Gn, are cal- 
culated by subtracting the mean value and normalized by the width of the distribution. As known from 
this figure the distribution of normalized growth rates changes its functional form for different values 
of N. The probability density functions are converging clearly to the theoretical function, p(G, 1.5,0). 
It should be noted that this convergence is rather slow with respect to N and we can approximate 
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Fig. 4 Width of growth rate fluctuation of Japanese business firms as a function of number of workers, N, 
in log-log plot. T he d ashed line shows a theoretical curve given by the inverse power law, jy -0 ' 046 , which is 
derived from Eq. ( 10 1 with a — 1.05. The inserted figure shows cumulative distribution function of annual sales 
of about 1 million J apanese firms in 2005 with the dotted line showing a power law distribution with the same 
exponent a — 1.05. 



the transient functional form typically by a double-exponential (Laplace) distribution as we can find 
for the case of N = 10 2 . This slow convergence can be a reason for observation of double-exponential 
distributions in various observational reports [1,6,8,15,17,22,23,24,32 . 

III. The case of 2 < a: Similar estimation for denominator of Eq. ^ is valid and the ordinary central 
limit theorem can be applied to the numerator of Eq. (|6|) as the variances for {xj(t)} are finite. The 



expression of Eq. ( 11 ) is also valid in this case, however, the parameters are limited to a = 2 and f3 = 0, 
namely, the limit distribution of normalized growth rate is always p(G; 2, 0), which is the well-known 
Gaussian distribution with no long tail. 

It is interesting to consider a special situation that the growth rates {gj(t)} are distributed sym- 
metrically around 1. Then, we can derive a = 1 from the basic relation < gj(t) >= 1, and /3 in Eq. 



(12 1 is by symmetry. In such a case we can expect that the limit distribution of normalized growth 



rate converges to p(G; 1, 0) = 1/ |2tt(1 + G 2 )} from Eq. (111. 

Results for all these cases are summarized in the third column of Table 1. The limit distribution 
of growth rate is determined by the value of moment function for the growth rates of subunits. The 
important point is that the ordinary central limit theorem can be applied only in the limited cases 
of relatively small growth rates, the mean value of growth rate is less than 1 and the second moment 
of growth rate is less or equal to 1. For systems in the real world it is expected that the systems are 
nearly in the statistically steady state and the mean values of growth rates of subunits may take a 
value around 1. Then, as known from Table 1, the limit distributions of growth rates belong to either 
power laws or non-universal functional forms. 



4 An application to business firm activities 



Now, we apply these theoretical results for data analysis. Among various types of dynamical complex 
systems in the real world, business firms are attracting the attention of scientists because there are 
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Growth rate 

Fig. 5 Growth rate probability density function of large business firms with more than 300 workers (black 
triangles) in semi-log plot. The red d otte d line shows the theoretical function of an asymmetric stable distribu- 
tion given by p(G; 1.05, 0.45) in Eq. (Jill) which is translated and re-scaled to fit the peak value and the width. 
The inserted figure shows the cumulative distribution function of the normalized positive growth rates (black 
line) in log-log plot to confirm the fitness with a theoretical stable distribution (red dashed line) and with a 
power law with tail exponent a = 1.05 (black dotted line). 



precise observation data in the form of financial reports pi[6ll MlTllT2lll"Tll2 1 II 3 3 j . In order to check 
the validity of our theory, we analyze comprehensive business firm data of Japan provided by the 
governmental research institute, RIETI (Research Institute for Economy, Trade and Industry). The 
data contains financial reports of 961,318 business firms, which practically covers all active firms in 
Japan in 2004 and 2005. It is already confirmed that basic quantities of these business firms, such as 
annual sales, incomes and number of employees, follow power law distributions |13j . 

There are several quantities that characterize the size of a business firm, such as assets, number of 
employees, sales, and incomes. Among these quantities, we focus on sales because this quantity reflects 
the present scale of activity of a firm most directly. Also, we simply assume that the whole activity of a 
business firm is given by the sum of the activities of individual employees; namely, we regard X{t; N) 
in Eq. Q as the annual sales of a business firm with N employees in the i-th year. 

Neglecting both the additive term in Eq. Q and the change of number of employees in a year, we 
estimate the growth rate G(t; N) by the ratio of t + 1-th year's annual sale over the £-th year's sale of 
a business firm with TV employees. It is already confirmed from the data that the autocorrelation of 
growth rate averaged over all business firms is very close to implying that the growth of a business 
firm can be roughly viewed as an independent random growth process approximated by Eq. Q with 
negligibly small external force term. 

Categorizing the business firms by the number of employees, N, we can measure the width of growth 
rate distribution for each category. Fig. [3] shows ./V-dependence of the width of growth rate distribution 
in log-log scale. It is confirmed that Gibrat's assumption of size-independence does not hold in this 
case, and the width of growth rate decays clearly for large N. Here, the theoretical line is given by a 
power law, jV~ - 046 , which is derived from Eq. (10) in the case of a — 1.05. In the inserted figure the 
cumulative distribution of annual sales of all firms, corresponding to a superposition of distribution 
of X(t; N) for all N, is plotted in log-log scale. We can confirm that the tail of the distribution is 
approximated by a power law with the exponent a = 1.05 as expected. 

A theoretical limit distribution of growth rate, p(G; 1.05, 0.45) in Eq. (11 1, is plotted together with 
that for real data estimated for N larger than 300 in Fig. [5] The inserted figure shows the cumulative 
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distribution of the normalized growth rate for positive G in log-log scale to confirm the functional form 
of the fat-tail. The growth rate distribution is asymmetric in this case and the whole functional form is 
well approximated by the theoretical curve of the stable distribution. This may be the first real-world 
example of application of an asymmetric stable distribution with a fractional value of characteristic 
exponent fitted in the whole range since the birth of mathematical theory in the 1930's. 



5 Discussions 



In this paper we introduced a theoretically solvable model of sum of randomly growing independent 
subunits. As summarized in Table 1 we found generalized central limit theorems applicable for a 
composite of randomly growing subunits, in which we can find various types in both width-size relations 
and the limit distributions of growth rates. As an example of a real world application, we analyzed a 
huge database of business firm growth rates of Japan and confirmed consistency with the theory. 

Application study of this theory for growth rate distribution of complex systems is highly en- 
couraged. It is expected that the shrinking of growth rate width and the functional form of limit 
distributions can be directly compared with real data of various phenomena to check the universality 
of this aggregated system of randomly growing subunits. As already commented for Fig. [T] the nor- 
malized growth rate distribution looks quite similar to a double-exponential (Laplace) distribution in 
the case of intermediate numbers of subunits. This type of finite size effect should be treated carefully 
in real-world data analysis. There is the possibility that the varieties of empirically known properties 
of the growth rate of complex systems introduced in the beginning of this paper can be understood 
using our approach as a frame work. 

Real-world systems may not be in a steady-state, so it is important for future work to investigate 
transient behaviors of this independent subunit system before it reaches the statistically steady state. 
Extension of this novel renormalization view of growth rates to cases of interacting subunits may also 
be an attractive new research topic. It is expected that a variation of generalized central limit theorem 
for the growth rates might be found for the wider category of growing complex systems like the case 
of non-standard statistical physics for long-range interaction systems |28j . 

Acknowledgements The authors appreciate RIETI for providing the business firm data. This work was 
partly supported by Research Foundations of the Japan Society for the Promotion of Science, Project Number 
22656025 (M. T.) and Grant-in-Aid for JSPS Fellows no. 219685 (H. W.). 



Appendix A: An introduction to random multiplicative process 

As random multiplicative process is not widely known, here we introduce a simple exactly solvable case of 
random multiplicative process and show intuitively how the process realizes a power law distribution in the 
statistically steady state. Also, a continuum limit version of this multiplicative process is discussed. 
We consider a positive random variable, x(t), which follows the following stochastic equation, 

x(t + At) = g(t)x(t) + 1, (A. 5) 

where g(t) is a stochastic noise term which takes either a positive constant, g, or with probability 1/2, 
respectively. Starting with the initial condition, x(0) = 1 , these time evolutions are given as 

*w = {? +1 K: I'M .«"•>- {P 

The general solution at time step r is obtained as 



1 (prob. 1/4) 
(prob. 1/4) . (A.5) 

(prob. 1/2) 



f-l)/(g-l) (prob. l/2 T ) 
r 1 -1) 1(9- 1) (Prob. 1/2 T ) 



x(rAt) 



-l)/(g-l) (prob. l/2 k ) 



(A.5) 



1 



(prob. 



1/2) 
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where k is an integer from 1 to r. From this solution we can calculate the cumulative distribution of x(t) in 
the limit of t — > oo as 

P(> x) = 2{1 + (g- l)x}"^M , (A.5) 

where P(> x) denotes the probability that x(oo) takes a value larger than or equal to x. In the case that g > 1 
we have the asymptotic power law distribution for very large value of x, 

P(>x)cxx- a , (A.5) 
where the exponent, a — log(2)/ log(g), fulfills Eq. Q in the range, < a < oo, namely, 

< g(t) a >= \ ■ g a + \ ■ o = i (A.5) 

With this special example we confirm the validity of Eq. |TJ to Eq. j3l. Note that in this example, the stationary 
condition, < log(<7(t)) >< 0, is automatically satisfied because the condition that g(t) = with probability 1/2 
gives the value, < \og(g(t)) >= — oo. 

The key point of realizing the power law in this multiplicative random process is understood intuitively by 
neglecting the additive term. The probability of repeating git) = g for k time steps is given as p(k) = l/2 fc = 
e -feio g (2) an£ j correspond ing v alue of x(t) is approximated as x(t) ~ g k = e fclog ' 9 ', then by deleting k from 
these relations we have Eq. (A.5 1. Namely, successive exponential growth with an exponential distribution of 
duration time gives the power law distribution. 

This type of derivation of power law can be generalized in the following way. Let us consider the following 
general form of random multiplicative process, 

x(t + At) = g(t)x(t) + f(t) (A.5) 

where g(t) and f(t) are inde pend ent random variables taking positive values, and we assume the situation that 
\ogg(t) fluctuates aroun d 130 ■ Taking logarithm for both sides and introducing variables y(t) = log (x(t)) 
and r(t) = logp(t), Eq. (A.5 I can be transformed as 

y(t + At) = log {g(t)x(t) + f(t)} = y(t) + r(t) + + • • • (A.5) 

g{t)x(t) 

Neglecting the terms including fit) as higher order terms, time evolution of the probability density of y(t), 
p(y,t), is approximated by a Fokker-Plank equation. 

p(y,t + At)n f oj(r)piy~r,t)dr^p(y,t)-<r>^p^ + ^^^^ + ---, (A.5) 

where uj(r) denotes the probability density of r. Assuming existence of a statistically steady state we have the 
following exponential distribution. 

2<r> 

p(y) oc exp<r' J > (A.5) 

In the situation, < r >=< log(g) >< 0, which is equivalent to the condition of existence of st eady state for 
random multiplicative process [T6], Eq. (A.5 1 is shown to be equivalent to the power law, Eq. (A.5 1, with its 
exponent given as 

a w t. — (A.5) 

This value is derived from the exact relation, Eq. (J3j) . By expanding the left hand side of the equation, 
< g(t) a >= 1, in terms a as follows, and by equating the second and third terms in right hand side as an 
approximation . 

< g(t) > = < e " > 

= 1 + a < lo S ( 9 (t)) > W <{1 ° S( f ))2}> + • - - (A.5) 



The key equation for determining the value of exponent, E q. j3t , can be derived roughly by the following way. 
Neglecting the additive term in the right hand side of Eq. ( |Ar5| and taking an average over realizations after 
taking the s-th power of both sides, we have the following relation. 

< x{t + At) s >^< g(t) s >< x(t) s >, (A.5) 

In the case that < g(t) s >> 1 it is clear that < x(t) 3 > diverges in the limit of t — > oo. On the other hand in 
the case that < g(t) s >< 1 the value of < x(t) s > is always finite. Therefore, we have the following relations, 

<*(~r>=fe4d. (A - 5) 
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Fig. B-l A schematic graph of the moment function. 



whe re a s atisfies Eq. Q3b . The property of Eq. (A. 5 1 is one of ty pical characteristics of the power law distribution, 
Eq. (A.5l. So, we can nnd that the power law exponent Eq. (A.5) is consistent with Eq. |3|. 

A rigorous mathematical derivation of this relation was done by Kesten in 1973 in more general form of 
real- valued matrix considering also the case that the distribution of the additive noise follows a power law [16] . 
In his proof the value of a is limited in the range of < a < 2 as he applies the theory of stable distribution, 
however, our numerical analysis and the above intuitive theoretical analysis suggests that the value of a can 
be extended to the whole range < a. 

It should be noted that the existence of the add itive term in Eq. ([I]) or Eq. (A.5 1 is essential to realize the 
statistically steady state. As known from Eq. (A.5 1 the stochastic process is a random walk with a negative 
trend in view of \og(x(t)), therefore, without any additive noise term the random walker tends to shrink to 
x{t) = 0, which is not the power law steady state. Even without the additive term (f(t) = 0) the same power 
law stea dy st ate can also be realized by introducing a repulsive boundary condition such as requiring x(t) > 1 
for Eq. (A.5 1 by adding a rule t hat x (t + At) = 1 when x(t) < 1. 

Dividing both sides of Eq. (A.5 I by At and considering the continuum limit of At — > 0, we have the 
following Langevin equation with a time-dependent viscosity. 



X(t) = -n(t)X(t)+F(t), 



where 



X(t) = x(t), n(t) 



lim 

At-ta 



i - g(t) 

At ' 



F{t) 



lim m. 

At-t-oo At 



(A.5) 



(A.5) 



As known from this equation the case of g(t) > 1 corresponds to a negative value of viscosity, fj,(t) < 0. In the 
case of a colloidal particle's diffusion in water such a negative viscosity cannot be realized, however, in the case 
of voltage fluctuation in an electric circuit, which is approximated also by a Langevin equation, we can consider 
a negative viscosity state by introducing an amplifier in the circuit. Namely, the value of fi(t) corresponds to 
resistivity in the electric circuit and in the situation that fluctuation of voltage is amplified as a whole, and 
the effective resistivity takes a negative value. By introducing an electric circ uit in which an amplifier works 
at random timing we have a physical situation w hich is described by Eq. (A.5 1 and power law distributions of 
voltage fluctuation are confirmed experimentally |28| . 



Appendix B: Basic properties of the moment function 

As Eq. |3| is the key of determining the exponent of the power law of Eq. |2|, a, here, we summarize the 
basic properties of the moment function for growth rate of a subunit, M(s) =< gj(t) a >. This is a continuous 
function and it is concave with respect to s for any distribution of gj (t) because the second derivative of this 
function is always positive, M(s)" =< (log(£,-(t)) 2 gj(t) s > > 0. As M(0) = 1 is an identity, if M{ot) = 1 holds 
for a positive valu e of a, then we know that M(s) < 1 for < s < a and M(s) > 1 for a < s as schematically 
shown in Fig. |B-l| So, M(2) < 1 corresponds to 2 < a as shown in the first column of Table 1. In the situation 
that Eq. (|3J) holds with < a < 1, then we have M(l) =< gj(t) » 1, while in the situation, 1 < a, we have 
M(l) XL 

The stationary condition, < \og(g(t)) >< 0, means that the slope of at the origin, M(0)', is negative, so 
if this condition is not fulfilled M(s) > 1 for any positive s, implying that the stochastic process of Eq. (fljl is 
not stationary. On the other hand in the case that the probability of occurrence of gj(t) > 1 is 0, it is trivial 
that M(s) < 1 for any positive s. 



Appendix C: A brief review of the generalized central limit theorem 



The central limit theorem is one of the most powerful mathematical tools; however, it is too often used to 
approximate the sum of random variables, of the form, Y(N) = yi + yi + • ■ • + yw, by normal distributions. In 
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fact there are three required conditions on random variables for their sums to obey the central limit theorem 

M- 

1. All random variables follow the identical distribution. 

2. The variables are independent. 

3. The variance of the variables is finite. 

If one of these conditions is violated, then the central limit theorem does not apply. 

Violation of the first condition has recently been attracting attention as super-statistics, that is, superpo- 
sition of stochastic variables having different statistics 5 j . As an old example, a fat-tailed velocity distribution 
observed in randomly stirred granular particles can be explained by superposition of normal distributions 
having different variances due to clustering caused by inelastic collisions |34| . 

Giving a general discussion of correlated variables violating the second condition is rather difficult, as the 
details depend on the details of correlation. It is known that universal properties independent of the details of 
the system can be expected at the critical point of phase transition at which power law distributions and power 
law scaling relations play important roles |31| . Also, theoretical approaches based on the concept of "non- 
extensive entropy" can provide general solutions for strongly correlated systems having scale-free interactions 
such as charged particles [13] . However, considering activity of a business firm, for example, it seems very 
difficult to describe both internal and external interactions by a general mathematical formulation. 

Violation of the third condition, infinite variance, was intensively studied in the 1930's by the pioneering 
mathematician P. Levy as "stable distributions" |19| . Assuming that {yi, 1/2, • • • , Un} are independent identi- 
cally distributed random variables, he showed that the fluctuation width of the sum Y(N) increases proportional 

to JV 1//a in general where a is called the characteristic exponent which lies in the range, < a < 2. The limit 
distribution is defined for the normalized variable, z = {Y(N) — cn}N~ 1 /°" , where cjv is a term corresponding 
to the mean value. In the limit of N — > 00, the distribution of 2 becomes a stable distribution that has a 
power law tail p(z) oc z~ a ~ x for < a < 2. The limit distribution converges to the normal distribution when 
a — 2, according to the ordinary central limit theorem. This general result is called the gene rali zed central 



limit theorem. The general functional form of the stable distribution is given by Eqs. ( 11 1 and (|12[) [9j 
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